// ----------------------------------------------------------------------
//  Jack McCaffery
//  11-30-2014
// ----------------------------------------------------------------------

#include "topicModel/multinomialParameters.H"

using namespace Eigen;

#include <iostream>
using namespace std;


#define VECTOR_SC Matrix<SC, Dynamic, 1>


namespace
{

/// \brief Project point_ onto the hyperPlane defined by the
///	   normal vector and a point in the hyperplane
/// \param point_ The point to project. Also the return value.
/// \param normal_ A unit normal vector to the hyperplane
/// \param coplanarPt_ A point in the hyperplane
template <class SC>
VECTOR_SC projectOnto (VECTOR_SC const &point_, VECTOR_SC const normal_,
		       VECTOR_SC const coplanarPt_)
{
	return point_ - (point_ - coplanarPt_).dot (normal_) * normal_;
}

}


///////////////////////////////////////////////////////////////////////////
template <class SC>
MultinomialParameters<SC>::~MultinomialParameters ()
{

}

template <class SC>
MultinomialParameters<SC>::MultinomialParameters () :
	ParameterGroup<SC> (),
	sumConstraint (1)
{

}

template <class SC>
void MultinomialParameters<SC>::applyGradient (
	FixedParameters<SC> const &gradient_, double rate_)
{
	if (gradient_.getData ().rows () != this->data.rows ())
		qFatal ("Data rows and gradient rows mismatch.");

	if (gradient_.getData ().cols () != this->data.cols ())
		qFatal ("Data columns and gradient columns mismatch.");

	unsigned lastIndex = gradient_.getData ().cols () - 1;

	Matrix<SC, Dynamic, Dynamic> constrainedGradient;

	// Constrain the gradient to the surface first. Then project into
	// the bounded region
	constrainedGradient = (rate_ * gradient_.getData ());

	constrainedGradient.colwise () -=
		(constrainedGradient.col (lastIndex));

	constrainedGradient.col (lastIndex) =
		(-1.0) * constrainedGradient.rowwise ().sum ();

	this->data += constrainedGradient;

	// Just keep projecting onto smaller hyperplane boundaries until
	// every dimension lands inside the bounds
	while (projectOntoBoundaries () > 0);
}

template <class SC>
void MultinomialParameters<SC>::setSumValue (SC sum_)
{
	sumConstraint = sum_;
}

template <class SC>
SC MultinomialParameters<SC>::sumValue () const
{
	return sumConstraint;
}

template <class SC>
int MultinomialParameters<SC>::projectOntoBoundaries ()
{

	VECTOR_SC normal (this->data.cols ());
	VECTOR_SC pointInPlane (this->data.cols ());

	unsigned numberProjected = 0;

	for (unsigned row = 0; row < this->data.rows (); ++row)
	{
		// The dimensionality of the hyperplane to project on to
		double dimensionality = 0;

		bool needsToProject = false;

		for (unsigned col = 0; col < this->data.cols (); ++col)
		{
			if (this->data (row, col) < 0)
			{
				normal[col] = 0;
				pointInPlane[col] = 0;
				this->data (row, col) = 0;

				needsToProject = true;
			}
			else if (this->data (row, col) > 0)
			{
				normal[col] = 1;
				pointInPlane[col] = 1;

				dimensionality += 1;
			}
		}

		if (! needsToProject)
			continue;

		numberProjected += 1;

		normal *= sqrt (dimensionality) / dimensionality;

		pointInPlane *= sumConstraint / dimensionality;

		this->data.row (row) = projectOnto<SC> (
			this->data.row (row), normal, pointInPlane);
	}

	return numberProjected;
}


template class MultinomialParameters<bool>;
template class MultinomialParameters<short>;
template class MultinomialParameters<unsigned short>;
template class MultinomialParameters<int>;
template class MultinomialParameters<unsigned int>;
template class MultinomialParameters<float>;
template class MultinomialParameters<double>;
